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We propose and study a continuum model for the dynamics of amorphizable surfaces undergoing 
ion-beam sputtering (IBS) at intermediate energies and oblique incidence. After considering the 
current limitations of more standard descriptions in which a single evolution equation is posed for 
the surface height, we overcome (some of) them by explicitly formulating the dynamics of the species 
that transport along the surface, and by coupling it to that of the surface height proper. In this 
we follow recent proposals inspired by "hydrodynamic" descriptions of pattern formation in aeolian 
sand dunes and ion-sputtered systems. From this enlarged model, and by exploiting the time-scale 
separation among various dynamical processes in the system, we derive a single height equation 
in which coefficients can be related to experimental parameters. This equation generalizes those 
obtained by previous continuum models and is able to account for many experimental features of 
pattern formation by IBS at oblique incidence, such as the evolution of the irradiation-induced 
amorphous layer, transverse ripple motion with non-uniform velocity, ripple coarsening, onset of 
kinetic roughening and other. Additionally, the dynamics of the full two-field model is compared 
with that of the effective interface equation. 



PACS numbers: 79.20.Rf, 68.35. Ct, 81.16.Rf, 05.45.-a 
I. INTRODUCTION 



Materials nanostructuring by ion-beam sputtering 
(IBS) has received increased attention in recent years^^ 
due to the potential of this bottom-up procedure for ap- 
plications in Nanotechnology, and also due to the inter- 
esting issues it arises in the wider context of Pattern For- 
mation at submicrometer scales^ In these experiments, 
a target is irradiated by a collimated beam of energetic 
ions (typical energies being in the keV range) that im- 
pinge onto the former under a well defined angle of inci- 
dence. Although routinely employed since long for many 
diverse applications within Materials Science (material 
implantation, sample preparation, etc.) the capabilities 
of this technique for efficient nanopatterninghave been 
recognized only recently, see references in [l||2|Q. Thus, 
it induces self-organized regular ripple (at oblique ion 
incidence) or dot (at normal ion incidence, or arbitrary 
incidence onto rotating targets) nanopatterns over large 
areas (up to 1 cm 2 ) on metallic, semiconductor, and in- 
sulator surfaces after a few minutes of irradiation. In- 
terestingly, the main features of this pattern formation 
process seem to be largely independent of the type of ions 
(even those inducing reactive sputtering) and targets em- 
ployed, as long as the latter amorphize under irradiation 
(the case of metals falls outside this class, and will not 



be addressed here, see e.g. [HHH). 

During IBS of amorphous or semiconductor substrates 
(for which the subsurface layer is amorphized, as fre- 
quently observed, see e.g. HQ) incident ions loose their 
energy through random collisions in the target bulkJ 
Near-surface atoms may receive enough energy and mo- 
mentum to break their bonds with the surface. Some 
of them may be certainly eroded, but most of them will 
be redeposited elsewhere, as seen e.g. in Molecular Dy- 
namics (MD) simulations^ In addition to adatom and 
vacancy formation^ which increases surface diffusion 
currents ^ enhancement of material transport by vis- 
cous flow seems to occur within a thin surface layer, as 
experimentally verified i 12 i 13 In any case, the evolution 
of the topography and the appearance of ordered pat- 
terns results from the balance between the erosive and 
the relaxational mechanisms. Whereas erosion tends to 
destabilize the surface (as a result of the fact that val- 
leys are eroded faster than crests^), relaxational pro- 
cesses tend to reduce height differences. Although there 
exists a wide separation of time scales between the hop- 
ping diffusive events (which are of the order of picosec- 
onds) and the ion-impact events (for an ion flux of 10 15 
ions cm~ 2 s _1 each atom on a typical surface experiences 
an ion impact about once per second) , both mechanisms 
have been modelled using kinetic Monte Carlo (kMC) 
approaches. The difference between both scales seems 
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to be fundamental to correctly describe the evolution of 
the irradiated surface — in typical time scales of the or- 
der of seconds 2 ^ — and challenges description by numer- 
ical simulations. In order to reach these length scales, 
a natural procedure is to resort to continuum descrip- 
tions. Hence, building upon Sigmund's description of 
the (Gaussian) energy distribution for energy deposition 
from collision cascades within amorphous or amorphiz- 
able targetsjii^, the seminal linear model of Bradley and 
Harper (BH)!^ and its non-linear extension o 16 i 17 ' 18 i 19 al- 
ready predict many of the experimentally important fea- 
tures, such as e.g. ripple formation and orientation as a 
function of incidence angle and dependencies of the ripple 
wavelength with temperature and flux. Moreover, they 
agree in many aspects with alternative models, such as 
kMC studies, see recent discussions in Q and Q. 

In all these continuum models a single evolution equa- 
tion is formulated for the surface height field, h(r, t), and 
contributions to such an equation are elucidated from the 
various relaxation mechanisms influencing surface topog- 
raphy. We will collectively refer to these as one- or single- 
field models. Nevertheless, they also present limitations, 
that we can group into several categories: 

(i) Inaccuracies of the energy distribution: The fact 
that there are known deviations from Sigmund's Gaus- 
sian distribution, most conspicuously at grazing angles 
of incidence^ may account for the incorrect order of 
magnitude of the roughening rate as estimated by these 
models, or their incorrect prediction^! for the direction 
of transverse ripple motion. 

(ii) Restricted number of mechanisms: Continuum 
models necessarily neglect physical mechanisms which 
may turn out to be important to the system behavior. 
This fact may be related with the unsatisfactory descrip- 
tion by one-field models of the ripple wavelength depen- 
dence with energy, phenomena such as pattern wave- 
length coarsening or, for the case of normal ion incidence, 
their lack of account for in-plane ordering, or the high pa- 
rameter sensitivity for dot formation. 

(Hi) Formal consistency: Under some circumstances, 
the very formal consistency of the one-field models can be 
questioned. For instance, due to the ad-hoc nature of the 
way in which competing physical effects (such as physical 
sputtering and surface diffusion) are merely added in the 
height equation of motion. Or due to the existence of 
cancellation modes of a varying natur o 19 ' 22 ! 23 ' 24 in the 
non-linear equations, or to physically unstable values of 
the effective surface diffusion coefficients for intermediate 
incidence angles J£ 

(iv) Non-linear features: Finally, the explanation for 
some of the experimental properties that remain insuffi- 
ciently accounted for by previous continuum models may 
require improvements on our understanding of non-linear 
effects (and thus, affect any further continuum descrip- 
tions). Some of these may include the direction of trans- 
verse ripple motion, the spread in the measured values of 
roughness exponents when there is kinetic roughening, 
and the value (as a function of physical parameters) of 



the saturated ripple or dot amplitude. 

Due to the insufficiencies of the current continuum de- 
scriptions of pattern formation by IBS, we conclude on 
the need for improved continuum models that (a) intro- 
duce increased number and/or type of relaxation mech- 
anisms in a natural way, that in particular allow assess- 
ment of the interplay between transport and morphol- 
ogy; (b) improve upon consistency issues (cancellation 
modes, etc.) of previous approaches; (c) can be adapted 
to modifications in the distribution of energy deposition; 
(d) can account for the phenomenology of nanopattern- 
ing by IBS within an unified framework, and (e) general- 
ize previous linear and non- linear models, incorporating 
their successes and improving upon their shortcomings. 

Trying to reach a balance between complexity and 
completeness in the physical description, in p5ll26ll27ll28j 
continuum models have been considered that are sim- 
pler than a full hydrodynamic description but still pro- 
vide an explicit coupling between the surface topography 
and the evolution of the relevant diffusive fields. Follow- 
ing the philosophy behind the so-called "hydrodynamic" 
approach to aeolian sand duneS ) 29 i 30 ' 31 in order to de- 
scribe the temporal evolution of the topography, two cou- 
pled fields are considered, namely, the density of mobile 
species being transported at the surface and the local 
height of the static target. Although naturally there are 
important differences between IBS nanopatterns and rip- 
ples on aeolian sand dunes (e.g. in IBS the size of the 
structures is roughly seven orders of magnitude smaller, 
the total mass is not conserved due to sputtering, and 
the nature of the morphological instability resides in the 
erosive process, rather than in the transport processes, 
as a difference with wind of water induced patterns on 
granular systems) both of them share global features that 
suggest modeling along similar lines. 

In this paper we study in detail this two-field approach 
to IBS, expanding previous results obtained in Refs. (27l . 
[28j . and focusing on the most generic case of arbitrary 
(oblique) angle of incidence that is pertinent to ripple 
formation. We will assess the extent to which two-field 
models can contribute to the improvement of continuum 
description of IBS as described in points (a) through (e) 
above, and can be seen as a continuum reformulation 
of thin film surface dynamics that goes even beyond the 
specific instance of IBS. Our aim here is also to clarify the 
influence of different experimental parameters, such as 
temperature or ion flux energy, in order to stimulate new 
controlled experiments. We derive an improved interface 
equation and relate the parameters appearing in it to 
experimental conditions and features. In a companion 
paper; 3 - 2 - that will be henceforth denoted as paper II, we 
explore the implications of our two-field model for the 
cases of normal ion incidence and rotating targets, that 
are of interest e.g. for the production of quantum dots by 
this experimental technique. 

This paper is organized as follows. In the next sec- 
tion the basic ideas of the coupled two-field model are 
discussed. In section IIIII its planar solution is obtained 
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and a linear stability analysis is performed. Section HV1 is 
devoted to obtaining a single effective evolution equation 
to describe the surface height of the bombarded surface, 
by means of a multiple scale analysis. In order to check 
the hypothesis made in the derivation of that effective 
equation, in Sec. [V] the dynamics of this equation will 
be compared with that of the original two-field model in 
the ID case. Following this, we will study the two di- 
mensional interface equation in Sec. IVIi and consider its 
relation to experiments. To end, we provide our main 
conclusions in Sec. IVIII In the appendices we collect de- 
tails of several analytical calculations. 



II. MODEL 

For the model formulation, a key experimental fact 
for amorphous and amorphizablc targets is the forma- 
tion through irradiation of a thin amorphous layer at the 
target surface, see references in 013 . As done in Refs. 
}25ll26ll27ll28| , the main model assumption is that the sur- 
face dynamics can be completely described through the 
time evolution of two fields: the height h(r, t) of the static 
substrate at time t and point r = (a;, y) on a reference 
plane that coincides with the uneroded flat surface, and 
the thickness R(r, t) of the (thin) surface layer of mobile 
species. This thickness can be related with the density of, 
say, mobile adatoms through their atomic volume. Note 
that for the energies we arc considering in the order of 
1 keV, we can take adatoms as the dominant diffusing 
species, although e.g. for energies below 1 keV, advacan- 
cies may dominate surface transport effects^ this should 
reflect in the values of the diffusion constant to be intro- 
duced below. 

Dynamics of the two fields are coupled, and read 



dtR 
dth 



(i - ^)r e 



Tad - V • J, 



(1) 

(2) 



where the x axis is chosen as the projection of the beam 
direction onto the xy plane. In H])-©, T ex and T a d, 
which depend on the geometry of R and h, are, re- 
spectively, the rate at which material is dislodged from 
the immobile target due to irradiation (locally decreas- 
ing the value of h), and the rate at which mobile ma- 
terial incorporates back into the immobile bulk (locally 
increasing the value of h). Therefore, in opposition to 
the excavation mechanism which is responsible for the 
overall decrease of h, there exists a process of incorpora- 
tion back to the bulk analogous of a local condensation 
of mobile species. Nevertheless, we will not consider a 
spontaneous rate of "evaporation" that is independent of 
the ion beam, so that we are neglecting surface tension- 
mediated evaporation/condensation effect a 33 ' 34 (equiva- 
lently, we are assuming that the pressure in the vapor 
phase is negligible). The excavated material may be ei- 
ther sputtered away, or added back to the mobile thick- 
ness R with an efficiency (1 — 0) = cf>. Therefore, the 



fraction of the eroded atoms which are finally sputtered 
away is represented by <f> so that, for 0^1, local redepo- 
sition is partially allowed^ For = 1 all eroded atoms 
are sputtered away, while in the = case the sputter- 
ing yield is zero. In the last case the effect of the ion 
beam is limited to providing material for surface trans- 
port, and there is no average motion of the interface. 
We will refer to the latter two cases as zero-redeposition 
and complete redeposition limits, respectively. They will 
constitute useful limiting cases below. 

The system ©-© was put forward in Refs. [HHfil, in 
which a linear stability analysis was performed. However, 
one of the limitations of the choices made in these works 
for T ex and T ac i is that surface diffusion vanishes in the 
absence of redeposition = 1, making the ensuing model 
ill defined (due to a large wave- vector instability). These 
limitations were overcome in Refs. [27.28], in which more 
physical mechanisms of erosion and addition are consid- 
ered. 

The third term on the right hand side of Eq. ([1]) de- 
scribes transport of mobile material onto the surface in 
the form of a continuity equation. In contrast to [2511261 ] . 
where terms representing Eiiich-Schwocbcl barrier effects 
(relevant to IBS of metals 1 -) are incorporated into the dif- 
fusive current J, these are not considered in (27H2H ]. With 
the aim of studying amorphous or semiconducting targets 
we will follow the latter option. Here we simply consider 
a diffusive term for mass transport onto the surface that, 
in the case of isotropic materials, is given by J = —DVR, 
where D may be a temperature dependent constant (see 
below) . 

Likewise, we will neglect momentum transfer in the di- 
rection of the projection of the beam of ions to adatoms, 
as this is expected to be non-negligible only at higher en- 
ergies (say, above 10 3 keV, see e.g. [36[ and a discussion 
in0). 



A. Excavation 

We next need expressions for the excavation and addi- 
tion rates. As studied in previous theoretical single-field 
studicS ) 15 i 16 ' 18 the rate at which material is sputtered 
from the bulk depends on experimental conditions such 
as the angle of incidence, 9, substrate and ion species, ion 
flux, <E>, average ion energy, E, temperature, T, and other. 
In these works, such dependencies were studied through 
an assumption on the shape of the spatial distribution 
for energy deposition, mostly Sigmund's Gaussian dis- 
tribution. However, there are cases in which systematic 
deviations from the Gaussian shape occur (see e.g. [37| 
for the occurrence of exponential decay combined with 
null energy deposition along the ion track). As recently 
shown moreover^ the shape of this distribution may af- 
fect the very existence of a morphological instability and 
thus the formation of a pattern. At any rate, given the 
fact that for most ripple patterns the aspect ratio is small 
enough so that a small slope approximation is expected 



4 



to hold 2 - (except, possibly, for compound materials and 
predesigned substrates^) , to lowest non-linear order, the 
form of the excavation rate is expected to b o 27 ' 28 

T ex = oo[l + a lx d x h + V • (a^Vh) + d x V ■ {g£7h) 
+ V • («4VV 2 /i) + d x h V • (asV/i) + Vh ■ K V/i)] (3) 

independently of the assumed energy 
distribution i 15 ' 16 i 18 ' 37 Here, we will ignore the ef- 
fects of direct erosion (knock-on sputtering) which could 
be relevant under very shallow energy deposition condi- 
tions {e.g. at very grazing angles of incidence). Indeed, 
the local erosion velocity that follows from Sigmund's 
distribution has the shape given in ([3]), see |l8[ and 
Appendix Changes in the energy distribution are of 
course expected to modify the values of the parameters, 
but not the number and shape of the terms appearing 
in ([3]), that are a consequence of the loss of x <-» — x 
symmetry induced by the oblique beam. Note that 
reflection symmetry is not lost in the y direction, and 
that the x — y symmetry can be restored under different 
incidence conditions, such as normal incidence {6 = 0) 
and for rotating targets, see paper II. Thus, we have 
that in general on = &ia,g{cti Xl ai y ) are 2x2 diagonal 

otixx o- ixy 



matrices for i = 2, 3, 5, 6, while a± 



The parameter ao defines the excavation rate of a flat 
surface and is directly related to the sputtering yield of a 
flat surface, Yq, the ion flux, and the number of atoms 
per unit volume in the solid, n v , by ao = <$>Yo/n v . Since 
typical fluxes range from $ = 10 12 cm -2 s -1 to <!> = 10 17 
ions cm~ 2 s _1 , the number of atoms per unit volume 
for an atomic diameter of 0.4 nm is n v = 30 nm -3 , and 
typical yields for experiments with ion energies of some 
keV are of order unity, then ao ~ 10 -3 — 10 2 nm s _1 . 

While the detailed dependence of the remaining a.yi 
coefficients on the physical parameters can be rather non- 
trivial, the main physical content of Eq. (J3j) is relatively 
straightforward. Thus, e.g., as already shown by BH, the 
coefficients aix and ct2 y are positive (see Appendix IA"| at 
small angles of incidence, which implies faster excavation 
at surface minima than at surface maxima, which is the 
landmark of Sigmund's morphological instability. Simi- 
larly, the various terms in ([3|) imply geometrical depen- 
dencies of the excavation rate with surface morphology; 
say^ for small 9 one has a± x < so that the excava- 
tion rate is larger on a lee (d x h < 0) ripple slope than 
on a stoss (d x h > 0) ripple slope. However, we will see 
that, when coupled to surface transport, some of these 
dependencies can be modified with respect to the sim- 
plest expectations. Conspicuous geometrical dependen- 
cies of this sort appear through the coefficients of the 
a4 tensor. Within Sigmund's energy distribution, these 
are high order geometrical dependencies of the sputter- 
ing rate that in one-field equations reflect into terms with 
the shape of surface diffusion. However, the present for- 
mulation makes it transparent the extent to which such 
terms do not correspond to actual material transport on 



the surface. We will come back to this point later. 



B. Addition 

One-field models are basically complete once T ex is 
provided. However, in our case we still need to specify 
the addition rate T ac i- To this end, we have to take into 
account that surface diffusion is an independent physi- 
cal mechanism that can take place even in the absence 
of an ion beam. Of course it should be susceptible of 
enhancement by the presence of the latter due to the 
induced increase in the density of diffusing species, but 
within our framework we would like to have surface dif- 
fusion currents which are not necessarily proportional to 
the ion flux. To this end, we will allow for a non-zero 
thickness of mobile material R eq even in the absence of 
excavation (T ex = 0) or redeposition (<f> — 0), and write 
down a rate that favors addition in highly coordinated 
surface positions (minima) rather than at sites with low 
coordination (surface maxima). Thus, we write^ 7 - 



7o -R 



Reqfi 



I2,d 2 x h ~ l2y d 2 y h)} 



(4) 



that has a form that is reminiscent from the Gibbs- 
Thompson expression effect for surface relaxation via 
evaporation-condensation. 33 ' 34 In Eq. ffl, 70 is the mean 
nucleation rate for a flat surface, 7Q -1 representing the 
typical time between two nucleation events, typically in 
the range of picoseconds, and 72^, 72^ > describe the 
variation of the nucleation rate with the surface cur- 
vatures. In principle this paper focuses on amorphous 
or amorphizablc surfaces, for which 72^ = 721/ = 72 
although, for the sake of generality, we will consider 
the most general case of anisotropic nucleation rates 
{l2x 7^ l2y) as far as convenient. 

As we will see later, the thickness of the mobile mate- 
rial, R, is only slightly altered off its equilibrium value, 
R e „, so that the rate of addition previously considered in 
[28j is equivalent to (01 , at least sufficiently close to the 
instability threshold. We will see in the next section that 
(IH) indeed leads to proper surface diffusion effects, that 
will allow us to identify the phcnomenological parameters 
R eq , D and 72 with physical constants. 



III. PLANAR SOLUTION AND LINEAR 
STABILITY ANALYSIS 

The existence of a wide separation of time scales be- 
tween diffusive events and erosive events will allow us to 
simplify the study of model (U)-([l|. We can assume than 
the excavation rate, ao, is much smaller than any other 
velocity involved in the problem. Specifically, by con- 
sidering ao *C joReq, we can define a non-dimensional 
parameter e = ao/(7oi? e q) which will simplify the study 
of the system in the following sections. As noted above, 
typically ao « 10 -3 — 10 2 nm s _1 , while the frequency 
of hopping diffusive events, equivalent to 70, is of the 
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order— of 10 9 s" 1 . If wc consider that the thickness of 
the mobile layer in equilibrium is of the order of some 
atomic sizes, R eq w 1 nm, we get as an estimate for typi- 
cal values of e to be in the range e « 10~ 12 — 10~ 7 , larger 
values corresponding to higher fluxes and/or larger yield 
conditions. 



A. Planar solution 

In order to start the study of our model, we first con- 
sider the situation of a perfectly flat interface. In such 
a case, all the spatial derivatives of h(r,t) are zero, Eqs. 
(P) and j2]) becoming 



d t R p 
d t h p 



: -ej Q R eq 



- 70 (R p - 

-io(R p 



Req) j 
Req) j 



(5) 
(6) 



where wc have defined R p (t) and h p (t) as the planar solu- 
tion fields. Integrating Eq. ([5|) and assuming R(t = 0) = 



R eq we obtain R p , which reads 



R p (t) = R eq [1 + £0(1 - e- 7nt )] 



(7) 



for any value (not necessarily small) of e. In |(7J| we see 
that, after a short time (of the order of Jq 1 ), R p reaches 
a stationary value equal to R eq , plus a small modification 
of order e due to the redeposition of excavated material 
(such an extra term is absent in the zero redeposition, 
4> = 1, case). As indicated in Sec. Ill Bt even in the ab- 
sence of excavation (ao = e = 0) or redeposition (0 = 0), 
there still exists an intrinsic fraction of mobile material 
equal to R eq . 

Substituting (|7|) into j6|) and assuming that h(t = 0) = 
0, we obtain the evolution of the planar height of the 
bombarded surface, namely, 



h p (t) = eR eq [-0 7o t + 0(£ 



-7o* 



-v t, (8) 



where the last expression holds for times longer than 7^ 1 , 
for which the planar profile erodes with a constant veloc- 
ity vo = e(j>joR eq = 4>otQ. This expression gives a clear 
interpretation of the parameter 4> as the overall efficiency 
of the sputtering process. 



B. Linear stability analysis 

The next step is to perform a linear stability analysis in 
order to investigate whether a small perturbation of the 
planar solution is amplified or damped out in the course 
of time. We consider periodic perturbations of the form 

i?(x, t) = R p (t) + R £ exp(ik • r + u k t), (9) 
h(x, t) = h p (t) + hi exp(ik • r + w k t), (10) 

where k = (k x , k y ) is the wave vector of the perturbation 
and <x>k its dispersion relation. Substituting Eqs. @ and 



(fT0| into (fT|) and ([2]), and neglecting quadratic terms in 
Rq, h,Q, we obtain the following linear system of equations 



where 



7o + Dk 2 
-7o Wk 



ex 

rl 



r 



Rn 



_ 



= 0, (11) 



=€ 7o i?. 



':'</ 



E 



OL\ x ik x - ^2 ( a 2j + a 3j ik x ) k 



r L = - lvReq{l2xk 2 x + l2yk 2 y ). 



(12) 
(13) 



Non-trivial solutions only exist when the determinant of 
the coefficient matrix equals zero, which allows us to ob- 
tain the dispersion relation Wk as the solution of the fol- 
lowing complex second order equation 



+ uj k (a + ib) + (c + id) = 0, 



(14) 



where the coefficients a, b, c, and d are functions of pa- 
rameters and wave- vector components, and are given in 
Appendix [Bl Eq. Q14p leads to two branches in the dis- 
persion relation, corresponding to its two (complex) so- 
lutions, namely, 



Re(c^) = -| ± -L ([(a 2 - 6 2 - 4c) 2 + (2ab- Adf] 



21 1/2 



Im(cj 



V2 

a 2 -6 2 -4c) 1/2 , (15) 
[(a 2 -b 2 -4c) 2 + (2ab~4d) 2 ] 1/2 

1/2 ■ (16) 



2 2V2 



-a 2 + b 2 



4c 



Substituting Eqs. (|B1|) - (|B4[) for a, b, c, and d into Eqs. 
(TTSj) and ([T5)l . we obtain an analytical expression for the 
dispersion relation as a function of the model parameters. 
Thus, we can describe the linear evolution of a periodic 
perturbation to the planar solution, since the real part 
of cjk is related to the growth or decay of the pertur- 
bation amplitude, while the imaginary part describes its 
in-plane propagation. Since we are interested in the be- 
havior of the system for long distances, we will reduce 
our analysis of Wk to small wave vectors. In this limit, 
we get, to lowest order in k x and k y , 



R- c Hc) = ~7o, 

Re(w^) = e07 O i? eg (a 2x kl 



a 2v k 2 y 



ecj)ai x kl) 



(17) 
(18) 



Thus, the negative branch is unconditionally stable (per- 
turbations decay exponentially for any wave vector) and 
non-trivial dynamics (including the pattern formation 
process) are thus governed by the positive branch, which 
features a band of unstable modes (wave vectors), of 
small magnitude for small e values, for which pertur- 
bations can grow exponentially. The imaginary part of 
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the dispersion relation for k = |k| -C 1 is lm(u>\^) = 
-e4>j Re q ai x k x , or Im(w k ) = ~e(jry R eq ai x k x depend- 
ing of the branch, the sign of k x , and the sign of (0—1/2). 
The linear in-plane propagation of the perturbations is 
related to the imaginary part of the dispersion relation 
[see Eq. (f2"9")) below]. For the positive branch, positive 
modes and large redeposition (0 < 1/2), we have 



Im(u£) = -ecj)-f R eq ai x k x , 



(19) 



which indicates that the perturbations travel along the a; 
direction with a constant velocity equal to olq4>ci\ x . 

With respect to the time evolution of the amplitude 
of perturbations, the linear pattern features are provided 
by that mode which makes Re(wj^) a positive maximum. 
For, say, small angles of incidence, both a 2x and a 2y 
are positive^ so that (JT5J) is maximized for infinite wave 
vector components. A finite maximum is seen to occur 
once we take into account higher order corrections (in k) 
to (|18p . where stabilizing mechanisms compete with the 
erosion instability. Thus, 



Re(w^) = e(jry R eq (a 2x kl + a 2y k 2 - e§a lx k 2 x ) 



Re q Dk 2 (j 2x kl+j 2y k;) 



HoReq y~] \<j)auj - { — - <t>R 
L V 7o 



Uql2i I 0*2] 



k i kj . 



(20) 



where we have kept terms that are lower order than 
0(e 2 fc 4 ). Eq. p0|) has the same form as the correspond- 
ing expression in one-field theories but with modified co- 
efficients, see Appendix [Cl In general, the 0(k 2 ) terms 
are both of a purely erosive origin, being directly propor- 
tional to the curvature dependencies of the excavation 
rate (once we neglect the 0(e 2 ) contribution), that are 
available for several energy distribution functions! 15 ' 18 i 37 
Thus, in particular, our model respects the signs of these 
terms as obtained e.g. by BHJ^ Given their destabilizing 
nature, they are usually referred to as "negative" surface 
tension terms. The remaining 0(fc 4 ) terms in (|20|) are of 
an opposite stabilizing nature related to surface diffusion 
effects as justified below. 



1. Two-field description of surface diffusion 

In order to clarify the physical meaning of the 0(k 4 ) 
contributions in (f20|l . it is useful to consider different 
relaxation mechanisms that are known to lead to such 
type of terms. 

a. Thermal surface diffusion Let us study at this 
point the extreme limit of no erosion in the original model 
(H|)- This can be achieved by simply "turning off" 
the ion beam flux setting ao = 0, which in turn im- 
plies e = 0. Note that, physically, in this case we are 
left with a system in which variations in the substrate 
height h are only due to local detachment/addition and 



transport of the surface mobile species R, precisely as in 
Muffins' classic description of surface diffusion activated 
by temperature . 33 i 34 Mathematically, the dynamics of 
the ensuing system ([I])-© conserves the total amount of 
material and, moreover, dynamics are linear (note, non- 
linearities enter only through the rate T ex , that has been 
turned off). Thus, one can readily solve the full system in 
this case. To our purposes we are interested in the long 
wavelength limit, for which we can simply take the e — > 
limit in the results of the present Section. Up to order 
0(fc 4 ), and already restricting ourselves to the isotropic 
case 723; = 72y = 72, the result is 



Re(w+) = -R eq D l2 k\ 
Im(w+) = 0. 



(21) 
(22) 



Thus, the exact evolution equation for the surface height 
in this long-wavelength limit reads 

d t h = -R eq D l2 V% (23) 

to be compared with Mullins' result ; 33 ! 34 namely, 



d t h 



k B Tnl 



V 4 /t, 



(24) 



where D s is the surface diffusivity of mobile surface 
species, v s is their concentration, 7 is the surface free 
energy per area, ks is Boltzmann's constant, and n" 1 
is the atomic volume. From this we see that the corre- 
sponding contribution in (|20[) is a generalization of sur- 
face diffusion in which the surface free energy is taken 
to be anisotropic in the two substrate directions. More- 
over, with the use of dimensional arguments, we identify 
parameters in T ac i as 72 = 7/(fcsTn„), D = D s , and 
R eq = u s /n v , whereby T ae i becomes an implementation 
of Gibbs-Thompson formula.— In any case, we see that 
the applicability of the two-field approach goes beyond 
the specific case of erosion by IBS, and it can serve as an 
intuitive phenomcnological reformulation of other phe- 
nomena within Surface Science. 

b. Surface confined viscous flow It is also a classic 
result^ that viscous flow, when confined to a thin surface 
layer, leads to a contribution to the height evolution of a 
similar form to (l24l) 



dth 



d 3 j 



(25) 



where d is the thickness of the viscous layer and rj s is the 
viscosity. In the case of IBS erosion of silicon targets, the 
relevance of such type of relaxation mechanism has been 
pointed out.— Specifically, it is argued in [13] that the ion 
beam induces this type of flow in such a way that 1 jr\ s cx 
E$, where E is the average ion energy. Notice that, 
under this assumption, all 0(k A ) terms in l]20[) would 
become proportional to ion energy and flux. 

In general, one expects both effects, thermal surface 
diffusion and ion-induced surface viscous flow, to occur 
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simultaneously in IBS systems^ 1 - so that an equation like 
(f23f should account for the effects described by (|24p and 
(|23|) . A form to accommodate this fact is to assume on 
a phcnomenological basis that R eq and D include both 
thermal (i.e. beam independent) and beam dependent 
contributions. 



2. Features of the linear instability 

We now come back to the full IBS model (i.e., for 
generic ao 7^ 0). Note that there are up to three dif- 
ferent 0(k 4 ) terms [second and third line in Eq. ([20]) ]. 
Besides thermal surface-diffusion of the type discussed 
in Sec. IIIIB 11 the terms proportional to otaj [on the 
last line of Eq. (12011 ] originate in the high order depen- 
dence of the excavation rate T ex with the height deriva- 
tives, and correspond to the so-called "effective smooth- 
ing" terms in one-field models. 18 ' 42 As is clear from our 
present formulation, being independent of R eq and D, 
these terms do not originate in actual material trans- 
port on the surfaced In marked contrast, the remaining 
0{k A ) terms on Eq. ([20)1 do couple excavation (they are 
proportional to ct2i) to surface transport (being propor- 
tional to cither D or i? e(J 72i), a feature that is beyond 
one-field descriptions. In particular, they may become 
temperature-dependent through the latter parameters, 
which will have relevant implications below. Similarly 
to one-field models, "surface diffusion" like terms oppose 
the erosive instability and lead to selection of a typi- 
cal length-scale in terms of the wave vector which grows 
(linearly) fastest. From (|20|) . we can obtain the features 
(orientation and magnitude) of such mode providing the 
ripple structure. 

a. Ripple orientation Using the results in Appendix 
[Cl for the small physically relevant values of e, the rip- 
ple structure can only align along the x or the y direc- 
tions. Using (|20| . for isotropic thermal surface diffusion, 
72a; = l2y , the ripple pattern is oriented along the x direc- 
tion (with crests aligned in the y direction) if a.2 X > &2yi 
or in the y direction (with crests aligned in the x direc- 
tion) when a.2y > Oi2 X , or is a linear superposition of the 
two orientations when ct2 y = ct2 X , in which case one has a 
square-symmetric cell arrangement, rather than a proper 
ripple structure. These results for the ripple orientation 
generalize those of one-field models ; 15 ' 18 for which there 
is moreover abundant experimental confirmation, see ref- 
erences e.g. in @0|. When thermal surface diffusion is 
anisotropic, 72a; 7^ "f2y, the possibilities of alignment for 
the ripple pattern are again along the x axis, along the y 
axis, or simultaneously in both directions (corresponding 
to an array of rectangular cells) if ctfxTss/ = a 2yl2x- 

b. Ripple wavelength In the cases above, the leading 
contribution (in powers of e) of the wave vector at which 
the linear dispersion relation is maximized reads 



"x,y 



U(jryoa2x,y 

2Dj2x,y 



1 ao4>a2x,y 

2R eq D-/2x,y 



where the x (resp. y) subindex applies when the ripples 
align in the x (resp. y) direction. Recalling the order 
of magnitude of the model parameters as given in the 
previous section, we can substitute them into (J2U). As- 
suming further <X2x,y an d l2x, y to be of the same order 
of magnitude {e.g., assuming that the only relaxational 
mechanism is thermal surface diffusion and employing 
the relations given in Sec. IIIIB 1 a[) . we have a.2x,y = 0.18 
nm and 72 = 3.8 nm using data for Si(001) as in |44J] fo r 
T = 500°C, and D « 10 5 nm 2 s" 1 as measured in [4(J. 



We thus obtain k 



10" 4 - 10- 



R, 



-1/2 



nm 1 , where 



we have used values for ao = 10~ 3 — 10 2 nm s _1 as above 
and the thickness of the mobile surface species layer, R eq , 
must be given in nm. If we consider this thickness to be 
comparable to a few atomic diameters, R eq sa 1 nm, we 
finally obtain an estimate of the linear ripple wavelength 
l l = 2tt /k l . Thus, r w 10 — 10 4 nm, in agreement with 
the experimental orders of magnitude^ 

Subdominant contributions to the ripple wavelength 
are physically very informative of the interplay among 
the physical mechanisms present in the two-field model. 
Thus, for instance in the case of ripples along the x di- 
rection one gets to next order in e 



I 1 = 2 3 / 2 7T 



2 / DR eq j2x A x Ct 4a 



1/2 



V (j)a a 2 x 4> a 2 x 
where we have used the parameter combinations 



<j)D 
7o 



(f>Reql2i, i = x,y. 



(27) 



(28) 



(26) 



In view of the physical interpretation of the various pa- 
rameters entering Eq. ([27]) . we see that the argument of 
the square root in this expression is the sum of a temper- 
ature independent contribution (the term a^ xx /a2 X ) cor- 
responding to the ion-induced effective diffusion of Ref. 
[l8[ and terms which include both thermal and beam de- 
pendent contributions. Such a compound structure for 
the linear ripple wavelength coincides precisely with that 
employed by Umbach et al.— when showing the impor- 
tance of surface viscous flow in order to account for the 
experimental behavior of the ripple wavelength with flux 
and temperature. It also has the same shape as that pro- 
posed in 0, capturing in a phenomenological way vari- 
ous experimental observations. We again stress that for- 
mula (|27p is obtained within a linear approximation for 
which the ripple wavelength is a time independent quan- 
tity. Thus, if ripple coarsening takes place in a given 
experiment, the finally observed wavelength is expected 
to depart from the value given by (|27p . 

c. Velocity of transverse ripple motion A third pat- 
tern feature that we can extract analytically within linear 
approximation is the velocity for transverse ripple mo- 
tion. This is the velocity at which, say, a local minimum 
of the linear ripple structure travels across the substrate, 
corresponding to the phase velocity of a wave packet 4 s - 
Note that the imaginary part of the dispersion relation 
only depends on the x component of the wave-vector, so 
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that (linear) ripple motion takes place only in the x di- 
rection. In order to compute its velocity we simply have 
to take the ratio between the imaginary part of the linear 
dispersion relation and the wave- vector, evaluating at the 
maximum of the real part of u k . Thus, 



h. This will be seen to allow for an adiabatical (pertur- 
bative) elimination of R from the dynamics. 

We will use a frame of reference comoving with the 
planar solution [Eqs. ([7]) and ©] in order to investigate 
how the solution evolves around it. We write 



V 



Im(w fe ) 



ocix + „ n , (— <f><X3x + A x ai x ) . 
(' ) 

(29) 

In the case of one-field models, an analogous expres- 
sion is obtained, except for the new term proportional 
to A x , that appears here due to the coupling between 
erosion and transport. Note the importance of an anal- 
ogous term (that is proportional to the ion beam flux 
and whose final sign is opposed to that of the combined 
first and second summands in (|29p . sec e.g. Appendix 
IA"| in order to correctly account for the experimental 
direction of ripple motion, as stressed in [2l|. In this 
reference, thermal spikes were invoked in order to jus- 
tify such an extra contribution. In contrast, our present 
two-field formulation allows to obtain a similar correc- 
tion [e.g. in the analogous zero redeposition limit we get 
V e = a a lx - {4ir 2 a /{l e ) 2 )(a 3x + Reql2xOii x )}, with- 
out the need for mechanisms that differ from, say, linear 
collision cascades combined with surface transport. Nev- 
ertheless, as with the ripple wavelength, nonlinear effects 
can in general influence the observed velocity of lateral 
ripple motion, as seen in Sec. IVI A 21 



IV. NONLINEAR ANALYSIS AND EFFECTIVE 
INTERFACE EQUATION 

During the development of the morphological instabil- 
ity, a time is reached after which nonlinear terms can no 
longer be neglected and a nonlinear analysis is needed. 
Note that the band of unstable Fourier modes extends 
from k* = \/2k x down to k = 0, its size being con- 
trolled by the square root of the small parameter e, as 
seen in Eq. (|26|) . Moreover, the fastest growing mode 
k e is also proportional to e 1 / 2 . Thus, e -1 / 2 provides us 
with a characteristic length scale associated with the lin- 
ear instability and makes it natural to define slow spatial 
variables that are of order unity at the scale of the linear 
ripple wavelength, namely, X = e 1 ' 2 x and Y = e x ^ 2 y. 
Moreover, it is also possible to obtain a estimation of the 
time scales associated with the translation (the imagi- 
nary part of uj£) and growth (the real part of w£) of 
the linear instability. Thus, by substituting the value of 
k £ in ([1^)) and (|2U)) . the imaginary part scales as e 3 / 2 
and the real part as e 2 . Hence, analogously to the slow 
spatial variables, we can define two slow time variables, 
Ti = e z l 2 t and T 2 = e 2 t, associated with in-plane trans- 
lation and vertical growth, respectively. These natural 
variables will allow us to perform a multiple-scale analy- 
sis in order to obtain a closed equation for h using the fact 
that, near the instability threshold (namely, for small e 
values), R tends to its stationary value much faster than 



h 
R 



h p + h 
R p + R. 



(30) 
(31) 



The strategy consists in expanding h and R in powers 
of e 1 / 2 , substituting these expressions into Eqs. fT} and 
@, and solving to increasingly higher orders in e. Before 
doing that we will write Eqs. ([I]) and ([2]) in terms of the 
slow, X, Y, Ti, and Ti variables by means of the chain 
rule 



e 1/2 d x , 
-e^ 2 dy, 
d t = e^ 2 d Tl +e 2 d T2: 



d x 



to obtain 



e z ' 2 d Tl R + e 2 d T2 R = 4>T ex - T ad + eDV 2 R, 
e 3 / 2 d Tl h + € 2 d T Ji = -f ex + f ad , 



with 



Tad = 70 [R + C^egV ■ (^V/l)J , 

Pex = loReq {f? /2 a lx d x h 

+ e 2 [V • (o2Vft) + V£ ■ (oaVft)] 

+ e 5/2 [d x V ■ (a 3 Vh) + 8 x hV- (a 5 X7h)} 



(32) 
(33) 
(34) 

(35) 
(36) 



(37) 



J V • (o4VV 2 /i)| , 



(38) 



where we have used the value of the temporal derivatives 
of the planar solutions, RP and h p , given by ([5]) and ©, 
expressed all space derivatives in the slow variables, and 
defined 72 = diag(7 2a; , 72j/)- 

Expanding now R and h in powers of e 1 / 2 as 



(39) 
(40) 



n=0 



we seek to solve for the various orders R ni h n by substi- 
tuting the above expansions into ([33)) and (|36p. 

Note that from Eq. ([55]) and substituting the value of 
T a d we obtain 



R = -eR eq V 



(7 2 V/i) + — U 
— 70 v 



T ex + eDV 2 R 



e 3/2 d Tl R-e 2 d T2 R], (41) 
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which, together with the shape of T ex given by (|38| . in- 
dicates that, for any order n, the R n coefficient depends 
on terms of lower orders in the expansion of R and h. 
Thus, the terms obtained in the expansion of R can be 
substituted back into to get a closed equation for 
the evolution of h. While details of this procedure are 
given in Appendix [D] the resulting equation reads, in 
the original time and space variables, 

dth = d t h p + ^xdxh 



| Vi dth + {dihf + ^d'fdxh + & (d x h)(dfh) 

i—x 7 y 



i,j=x,y 



(42) 



where we have neglected height derivatives that are of 
sixth or higher orders, we have undone the transforma- 
tion to the frame comoving with the planar solution, and 
parameters are related to those of the original two-field 
model ©-© as 



Ix 



Oilx, 



v x = -a Q 4>a 2 x H -(j>(f>a\ x , v y = -a 4>a2 y , 

7o 



A 



(i) 



oco4>a 6l , d = -a (j)a 5 i, 



f2j = cto 
)Cij = DR eq ^ 2l + ao 



-<pOL3i + (j>Reqj2i aix 

V 70 



ba izj - I 4>Reql2i I a 2j 

7o 



\( 2 ) ( ^ D AO 



(43) 



Note that in (|43|) we have restored the expression of e in 
terms of physical parameters. As mentioned in Sec. lIII Al 
after a time of order ^y^ 1 the the profile erodes uniformly 
with velocity vq = —dth p = 4>cto- 

We have obtained a closed evolution equation for h 
from which R has been eliminated, and whose behav- 
ior is equivalent to that predicted by the full two-field 
model near the instability threshold. Note that, in par- 
ticular, the linear dispersion relation for (|42p coincides, 
within our long wavelength approximation, with that of 
the original model as given by (fl9| and ([20]) . Moreover, 
as in previous one-field descriptions, in Eq. (|42p there 
is not reflection symmetry in the x direction due to the 
oblique ion incidence. This symmetry is restored if the 
bombardment is perpendicular to the substrate, or else 
if the target is rotated simultaneously with irradiation, 
as described in (32j. Actually, Eq. (|4"2|) generalizes the 
anisotropic interface equation (|A1[) that is obtained by 
one-field theories^ by the appearance of additional non- 
linear terms (with coefficients A 2J ■ ). These, together with 
the modified dependence of parameters on physical con- 
stants, are the main effects of having explicitly described 
the dynamics of the diffusive field R onto the evolution 



of the profile, and will be seen below to be instrumental 
in order to provide an improved description of nanopat- 
tcrning by IBS. 



V. ID MODEL 

Eq. (|42p is a highly non-linear and anisotropic system 
whose full analysis is rather complex. Before analyzing 
it in detail, and in order to understand more directly the 
physical content of its various terms and parameter de- 
pendences on physical constants, we are going to study 
first a ID counterpart of the erosion model studied in pre- 
vious sections. We will thus consider that the x axis is 
the only relevant direction to describe the topography of 
the system. This simplification is very frequently done in 
models for sand ripples formation , 29 ! 30 ' 31 in which trans- 
lation invariance is assumed in the direction perpendic- 
ular to the wind. Note that such an approximation still 
respects the physically essential lack of reflection sym- 
metry in the x axis. Thus, by repeating the approach of 
the previous section in the case that there is no variation 
of the fields in the y direction, we obtain the following 
one-dimensional equation 

d t h = -no + ixdxh + uxd 2 c h + \^{d x hf + n x dlh 
+ ix (dxh)(d 2 x h) - Kxxdth - AS d 2 x (d x h) 2 , (44) 

where by an abuse of language we will employ similar 
symbols for parameters to those of the previous Section, 
and the relation of these with the coefficients of the cou- 
pled model are 



v = a o 0; 7a = -a (pai x ; v. 



-a (pa 2 x H #a l3 

7o 



4>Ot?,x + [ — 4>Reql2x J «la 



;Cx = -a (/)a 5x ; 



^ J. u 

(f>Reqj2x Oi2x 

7o 



JCxx = DR eq j2x + a 

= -a 4>a 6x ; A^ 2 J = -a (— - <t>Reql2x J a 6x 

V 7o / 



(45) 



Eq. (TJ3|) provides the generalization of the ID counter- 
part of Eq. (|Aip . through appearance of the additional 
\xx term. Actually, restricting ourselves to even terms 
in x derivatives (that is, for j x — ^x — Cx = 0), 
Eq. (|44[) becomes the mixed Kuramoto-Sivashinsky equa- 
tion (see [46| and references therein) that generalizes 
the Kuramoto-Sivashinsky (KS) equation . 47 ' 48 In general 
note that the coefficients (|45| directly reproduce those as- 
sociated with the x direction among the larger set of pa- 
rameters in (|43|) . Although one dimensional, Eq. (|44| is 
still a highly nonlinear equation with behaviors that may 
range from in-plane traveling periodic (ordered) struc- 
tures to chaotic (disordered) cell dynamics, as occurs 
with its \ { x 2 } = limit 
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A. Physical interpretation of parameters 

Before attempting to understand the interplay among 
the various terms in Eq. (|44f , it is worth giving consider- 
ation to each one of them individually. To this end, it is 
instructive to start by studying the two possible limiting 
cases for parameter <fi. 



1. Complete redeposition ((f) — 0) 

Equation (|44|) becomes strongly simplified when the 
erosive mechanism limits itself to transferring material 
from the immobile bulk to the mobile diffusive current, 
without sputtering proper, akin to the role of IBS for ion 
beam assisted deposition*^ In this case, the only non- 
zero coefficients in (|45|) are 

_ a Dai x (2) _ a Da 6x 

x ~~ „, ' 11 ~~ — ' 
7o 7o 

K. xx = DR eq j2x , (46) 

7o 

the interface equation reading merely 

d t h = n x d 3 x h - IC xx d 4 x h - d 2 x {d x h) 2 . (47) 

This equation has the conserved form expected from 
the fact that excavation is here limited to matter re- 
distribution. Actually, in the absence of the third order 
derivative term, Eq. (|4T|) in known as the conserved KPZ 
equation ) 52 i 53 i 54 relevant to conserved surface growth dy- 
namics such as in typical Molecular Beam Epitaxy sys- 
tems. Note that, although the surface diffusion coefficient 
IC XX of Eq. (|46p includes an erosive contribution that is 
of a destabilizing nature as long as excavation is favored 
at surface minima {a.2 X > 0), being proportional to uq 
this contribution is numerically smaller than the stabiliz- 
ing (thermal) contribution also present in JC XX . The only 
remaining nonlinearity in (|47p reflects (through ae x ) the 
non-linear dependence of the excavation rate with the lo- 
cal surface slope. Moreover, already this term genuinely 
couples erosion to transport, being also proportional to 
D. 



2. Zero redeposition (cf> — 1) 

This limit corresponds to the usual assumption in pre- 
vious one-field approaches. In this case gencrically Eq. 
(|44]l displays its full shape, with coefficients 

vq = a ; 7 X = -a ai x ; v x = ~a a 2x 

= -"0 ("3a; + Reqj2xOilx) £3 = -Oi Q a 5x ; 

ICxx = DR eq ^i2x + a>o (a ixx + R e qj2xOi2x) ; 

A^ 1 ' = -a a 6x ; X x 2 } = a R eq j 2 xaex- (48) 



Among coefficients in (|48[) . all but three of them (fl x , 

(2) 

K. xx , and X x i ) are directly as predicted by one-field mod- 
els, sec (|A2p . As for the three remaining coefficients, 
common to all three is that they correspond to conser- 
vative terms in the equation of motion. This allows to 
understand the contributions that they include in which 
transport (through dependence on R eq ) couples to an ero- 
sive dependence on a height derivative two orders lower. 
E.g. £l x is associated with a third order height deriva- 
tive and indeed features a direct erosive dependence in 
the 3rd. order coefficient a^ x . However, it also depends 
(through R eq ) on the first order erosive coefficient a\ x . 

Similarly for 1C XX and X xx . The surface diffusion co- 
efficient IC XX adds to these the expected contribution 
DR eq ^2x discussed in Sec. IIIIB 11 Moreover, note that 
the ion effective smoothing term with coefficient a^ xx , 
that reflects the dependence of the excavation rate with 
high (fourth) order surface derivatives, appears as a di- 
rect contribution to the surface diffusion coefficient. 

About the coefficient of the conserved KPZ term, note 
that for this <f> = 1 case its sign is opposite to that of X x ^ 
in |48|) . This leads to a cancellation mode and math- 
ematically invalidates Eq. (|4"4"|) description of the 
physical system. Indeed, neglecting the £ x nonlinear- 
ity that does not participate in the height saturation of 
the linear instability^ the remaining nonlinear contribu- 
tions read, in Fourier space, -(A<b 1} + k 2 \ x 2 })T[{d x h)% 
where T denotes Fourier transform. Due to the signs 
of the coefficients, there is a wave vector in the unsta- 
ble band (cancellation mode) for which the parenthesis 
in this equation vanishes, rendering the system nonlin- 
early unstable^ This undesirable feature actually also 
occurs in full 2D one-field models when generalized to 
sufficiently high orders. 19 i 23 ' 24 



3. Partial redeposition (0 < (j> < 1) 

Generically we expect partial redeposition to occur un- 
der usual experimental conditions for IBS nanopattern- 
ing. After the previous Section, we see that not only is 
redeposition a physical effect to include but also that it 
allows to regularize our mathematical description of the 
system. Indeed, using the parameter combination de- 
fined in {28}, we see that parameter conditions exist for 
small but non-zero values < < 1, for which A x > 
so that Ax and Xxx have the same sign and cancella- 
tion modes do not occur. The numerical values of <fi and 
A x also affect the remaining coefficients in (|45p. but are 
of a less critical nature. The only contribution that is 
privative of these partial redeposition conditions is the 
second term in the expression for v x , that, being posi- 
tive, is of a stabilizing nature and opposes the sputtering 
instability. A similar term can be found in the forma- 
tion of macroscopic ripples under the action of the wind 
when the number of sand grains is not conserved^ and 
reflects the geometrical fact that erosion tends to smooth 
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out inclined surface features. Nevertheless, such a term 
being higher order in powers of ao, we expect it to be 
numerically small in most practical cases within our IBS 
context. In general, the < <j> < 1 case interpolates be- 
tween the two extreme cases considered above, in that 
the dependence of coefficients (|45p on physical param- 
eters combine the features discussed in Sees. IV A II and 
IVAl 



Effective interface equation vs full two-field 
model 



In order to check the analytical approximations made 
in the derivation of the effective interface equation and 
compare its predictions on the dynamics to those of the 
full original two-field model, we have performed a nu- 
merical integration of the ID coupled set of Eqs. (HJ- ([2j) , 
and of the related single Eq. (|4"4"1) , using an Euler scheme 
for the time integration, and the improved spatial dis- 
cretization introduced by Lam and Shin 5 ' for the nonlin- 
ear terms. We have used periodic boundary conditions, 
lattice constant Sx = 1 and time step 5t = 0.01, checking 
that results do not differ significantly for smaller space 
and time steps. The standard system size of our simula- 
tion has been L = 2048. With the aim of comparing the 
evolution of the profile for the two-field and the effective 
equations, the same random initial height values were 
chosen, uniformly distributed between —0.05 and 0.05, 
and the corresponding parameters were related following 
(H5J). 

We show in Fig. Q] the evolution described by the ID 
two-field model of the height profile h and the thick- 
ness of the mobile material above h for certain values 
of the parameters. Since e = 3 • 10~ 3 is small, we see 
that R is indeed only slightly altered from its equilib- 
rium value (R eq = 1). Note how the morphological in- 
stability leads to formation of a periodic ripple pattern 
that, as expected, is not symmetric in the x direction. 
The thickness of the mobile surface layer correlates with 
the topography all along the dynamics, being smaller at 
steeper ripple sides. 

In Fig. [2] we compare the evolution of the profile for the 
ID two field coupled model, with that described by the 
effective height equation, Eq. (|44l) . where the coefficients 
of both systems are related by (|4"5"j) . We can see how, 
starting from the same flat random initial distribution 
for both systems, a periodic surface structure appears 
with a wavelength of about the maximum of the linear 
dispersion relation, and the amplitude of height varia- 
tions increases. For the examples considered in Fig. [2] 
the wavelength of the linear instability is given by ([27]) . 
yielding I = 98. At these short times, when the slopes 
are not too large so that nonlinear terms are not yet rel- 
evant, both profiles match quite accurately. Far from the 
linear instability threshold, the profiles become less sim- 
ilar. Since the space and time scales separation and the 
power expansion performed to obtain the effective equa- 
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FIG. 1: Height profiles h (dashed line) and thickness of the 
mobile material over h (solid line) at different times given by 
the ID two-field coupled model with <j> = 0.99, ao = 0.03, 
oti = — 1 Q2 = 30, az — cti — 1, a 5 = — 1, ae = —3, R eq = 
7o — 72 = 1, and D — 10. All units are arbitrary. 
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FIG. 2: (color online) Height profiles at different times given 
by the ID two-field model (black line) for parameters as in 
Fig. [Hand by the effective equation [Eq. I|44[l] (green line) with 
parameters as given by relation (|45|) . namely, vq = —3 ■ 10~ 4 , 



7* = 3 • 10" 4 , v x = -9 • 10" 3 , = 9 • 10~ 4 , tl x = -0.297 , 
£x = 3 • 10" 4 , K. xx = 1.0993, and A^i = 0.8901. All units are 
arbitrary. 



tion are only valid for small values of e, it is expected 
that, the smaller e is, the more similar the profiles be- 
come. However, if we reduce this parameter, the simu- 
lations are more time consuming since the characteristic 
space and time scales for pattern dynamics are inversely 
proportional to powers of e, as noted in Sec. IIVI In any 
case, for the values of e considered in our simulations, 
the effective equation captures the main features of the 
original two-field model, even in terms of the behavior 



of observables such as the global surface rms width or 
roughness W(t) or the ripple wavelength l(t) (defined as 
the mean lateral distance between two consecutive local 
minima), as seen in Figs. [3] and [51 respectively. 



C. Nonlinear dynamics for the effective equation 

Indeed, at later stages, nonlincaritics determine the 
evolution of the surface morphology. For the reasons 
mentioned above, we will explore this regime through 
the effective interface equation. Specifically, nonlinear 
effects induce coarsening of structures wherein the cells 
(ripples) grow in width and height, their number decreas- 
ing in both systems. For both cases coarsening is such 
that smaller cells are "eaten" by larger neighbors un- 
til reaching constant amplitude and wavelength values, 
while lateral moundlike order is still preserved for inter- 
mediate distances (more than ten times the lateral size of 
the cells). This behavior is very similar to that reported 
in Ref. [46| for the mixed KS equation equation that cor- 
responds to the j x = Q x = £ x = limit of ([33]); see also 
paper II. 

In Figs. [3] and [4] we show the time evolution of the sur- 
face roughness W(t) and ripple wavelength l(t), averaged 
over 18 random initial conditions. After a stage in which 
the amplitude of the linear instability and, therefore W, 
grow exponentially, a coarsening process begins (roughly 
at t ~ 3T0 5 ) for the ripple wavelength. Around t = 2-10 6 
this process stops and the wavelength and amplitude of 
the pattern reach stationary values. Specifically, the lat- 
eral pattern wavelength grows from its initial value cor- 
responding to the linear instability I = 98 until a sat- 
uration value, close to I — 121. At intermediate times 
this coarsening behavior can be described by an effec- 
tive power law I ~ t ' 12 , as suggested in Fig. [4j In the 
presence of coarsening, the dependence of the asymptotic 
values of the ripple amplitude and wavelength with sys- 
tem parameters differs from those of the linear instability 
regime. If one assumes^ that the odd-derivative terms 
in Eq. ([44]) do not contribute to such a coarsening pro- 
cess, approximate values can be obtained through com- 
parison with coarsening dynamics in the conserved KS 
equation^ Such estimates are more accurate in the nor- 
mal incidence case (paper II) ^ to which we refer the in- 
terested reader. Additional important features of these 
systems, which are not present in the equation studied in 
Ref. [4g|, are the asymmetry of the profile and the lat- 
eral movement of the pattern. As we have checked in our 
simulations, the asymmetry on the pattern depends only 
on the (advective) terms corresponding to the coefficients 
Q x and £ x of the effective equation. For negative values 
of fl x and/or positive values of £ x the cell structure tends 
to be leaning to the right. This can be observed in Figs.Q] 
and[2] where the right slopes of the cells are clearly larger 
than the left slopes. If fl x is positive and/or £ x is nega- 
tive, the pattern is leaning to the left. If both terms have 
the same sign, the orientation of the structure depends 
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FIG. 3: (color online) Temporal evolution of the global rough- 
ness, W(t) given by the two-field model (blue triangles) and 
by the effective interface equation (green circles) for the same 
coefficients as in Fig. [2] Error bars are smaller than the sym- 
bol sizes. All units are arbitrary. 
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FIG. 4: (color online) Temporal evolution of the lateral pat- 
tern wavelength, l(t), given by the two-field model (blue tri- 
angles) and by the effective interface equation (green circles) 
for the same coefficients as in Fig. [2] A few representative 
error bars are given that represent statistical dispersion. The 
dashed line corresponds to l(t) ~ t° . All units are arbitrary. 



on their relative magnitude. 

Considering lateral ripple motion, note first that the 
linear prediction for the velocity, Eq. ([29]) has the form 
V 1 = —j x + 4ir 2 £l x (l )~ 2 . The contribution due to j x is 
an uniform translation (along a direction on the x axis 
that is opposite to the sign of j x ) that can actually be 
cancelled out by an appropriate choice of reference frame. 
Thus, the only remaining terms which influence in-plane 
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displacement of the pattern are again fl x (linear) and £ x 
(non-linear). For the values we have considered for the 
remaining parameters, a positive sign of f2 z and/or £ x 
induces ripple motion towards positive x, while negative 
values of these parameters lead to lateral ripple motion 
in the opposite direction. In Fig. [5] we observe the lateral 
movement of the pattern as described by Eq. (|44p . for 
parameters as in Fig.O Simultaneously with erosion and 
mean height evolution towards larger negative values, the 
pattern is moving towards the left. Here, the movement 
is dominated by 7^ and Cl x , which induce motion towards 
the negative x values. 
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FIG. 5: Height profiles between t = 1.5 • 10 6 (top) and t = 
1.8 • 10 6 (bottom) evaluated at equally spaced intervals of 10 4 
time units, for the effective equation (|44|) . and parameters as 
in Fig. [2] All units are arbitrary. 



VI. FULL 2D EFFECTIVE INTERFACE 
EQUATION 

Equation (|42[) generalizes the one-dimensional equa- 
tion (|4"4"| to the case of fully anisotropic two-dimensional 
targets, in a way that is consistent with reflection sym- 
metry in the y direction, as expected from the ion inci- 
dence geometry. As mentioned earlier, Eq. (|42[) general- 
izes the one- field equation (|A1[) , through appearance here 
of the (anisotropic) conserved-KPZ type terms df(djh) 2 . 
In turn, Eq. (|A1|) already provided an anisotropic gener- 
alization (through the presence of odd derivatives in the 
x coordinate) of the two-dimensional KS equation. 16 i 22 
To the best of our knowledge, Eq. (|4"2")l is new and adds 
to the relatively small number— of (local) evolution equa- 
tions for fully anisotropic two-dimensional pattern form- 
ing systems, that are derived from constitutive laws. In 
the context of hydrodynamic models of ripple formation 
on aeolian sand dunes, an isotropic 2D equations, when 
available, are limited to conservative dynamics,— while 
in the cases of thin film surfaces nonlinearities that arise 
are of a different type i 59 ' 60 i 61 

Although the parameter space of Eq. (|4"2"|) is much 
larger than that of its one-dimensional counterpart (j44[) , 
the physical interpretation of the various terms and coef- 
ficients is completely analogous, corresponding to a natu- 
ral generalization of those appearing in the latter. Given 
that the main linear features of the two-dimensional 
equation were already discussed (and compared with typ- 
ical experimental data) in Sec. IIII B 21 we next consider 
numerical simulations of Eq. ([42]) that show the main 
morphological features of its full dynamics, that will be 
later compared with experimental results. Some pecu- 
liarities on the cancellation modes that may arise in Eq. 
(|42|) are considered analytically in a specific subsection. 



The results reported in this section allow us to con- 
clude that both the effective interface equation and the 
two field model, whose coefficients are related through 
(|45p . capture common features observed in experiments 
such as the coarsening process of the pattern wavelength, 
the short range lateral order and the non uniform lateral 
displacement of the structure. On the other hand, due to 
the fact that the scales associated with the experimental 
linear instability are very large (of the order of e -1 / 2 ), one 
needs very large simulations in order to compare with ex- 
periments. These are available to the effective equation, 
in which parameters can be rescalcd with the aim of ac- 
celerating the simulations. For these reasons, in going to 
the physical 2D case in the next section, we will limit our 
study to the 2D effective height equation (|4"2"j) . We will 
consider some illustrative examples of the ensuing sur- 
face dynamics that allow us to understand the richness 
of the behaviors that can be described by such a complex 
non-linear system. 



A. 2D dynamics: numerical results 

Far from a complete analysis of Eq (|4"2")) , we will limit 
ourselves in this section to a qualitative study of its main 
properties and how it successfully reproduces some ex- 
perimental features which are not included in previous 
continuum descriptions. 

Thus, we have performed a numerical integration of 
Eq. (|42|) using an scheme that generalizes that employed 
in the one-dimensional case, namely, an Eulcr updating 
rule with St = 10~ 3 for the time evolution, and the fi- 
nite difference prescription of [57] for the nonlinear terms, 
with 5x = 1. The standard size of our simulations was 
L x L = 256 x 256 with periodic boundary and random 
initial conditions. Wc consider a reference plane comov- 
ing with the eroded surface with a constant velocity — vo, 
thus the effective equation that we integrate is (|42"|) for 
vo = 0. 

The evolution of the height as described by Eq. 
is depicted in Figs. OH1 for different values of the coeffi- 
cients, with the x-axis oriented along the horizontal di- 



rection (see also supplementary videos) i 62 i 63 ' 64 In each 
figure three snapshots (top views and lateral cuts) are 
provided for a given parameter condition, with time in- 
creasing from panel (a) to panel (c). In all these exam- 
ples, and resembling experimental morphologies^ both 
the amplitude and the wavelength of the ripples grow 
with time, while the pattern disorders in heights for long 
lateral distances. The detailed shapes of the topogra- 
phies, however, are quite different depending on the val- 
ues of the parameters. We can obtain longitudinally dis- 
ordered ripples which are frequently interrupted along 
the direction of the crests, as in Fig. O or else ordered 
straight and wide ripples occur for different parameter 
conditions as in Fig. [7l An even more disordered pattern 
is depicted in Fig. [8j where the ripples group into do- 
mains of about three cells whose crests run along the x 
axis, as expected from the parameter values (note v x > 
in this example). 

Similarly to the one-dimensional case, before slopes are 
large enough to make non-linear terms non-negligible, 
the evolution of the morphology is governed by linear 
terms. This will allow us to separate the dynamics into 
two different regimes, linear and nonlinear, according to 
the type of terms that control the evolution. 



1. Linear regime 

As noted in Sec. IIV( the linear dispersion relation of 
Eq. (|42|) coincides with that of the original model de- 
scribed in subsection IIII Bl Thus, for isotropic thermal 
surface diffusion, the ripple crests are oriented along the 
y (x) axis if v x (v y ) is more negative than v y (u x ), thus re- 
producing the ripple orientation as predicted by the BH 
theory. Numerical integration within linear regime in- 
deed retrieves the dependence of the ripple orientation as 
a function of the values of v x and v y as shown in Figs. El 
El and[5] Furthermore, we have also checked in our simu- 
lations that the lateral wavelength of the pattern is given 
by the relation between the surface tension and diffusion 
terms. One way to do that is to measure the distance 
from the origin to the hrst maximum in the height auto- 
correlation function which is represented in the inset of 
Figs. EJa),[7Ja), andEJa). Since /Cy = 1 is considered for 
all these examples, we have l e = 27r/fc ? £ = 27r(— 2/z^) 1 / 2 . 

While even linear derivatives in Eq. (pf!?)) are respon- 
sible for amplification or attenuation of the ripple am- 
plitude, they do not induce lateral motion of pattern. 
Conversely, odd derivatives breaking the x — > — x sym- 
metries indeed induce in-plane lateral ripple motion. We 
have checked in our simulations that, as expected, the 
term corresponding to the coefficient "f x does not alter the 
shape of the morphology but merely produces a uniform 
movement along the x axis. As in the one-dimensional 
case, the direction of this movement is opposite to the 
sign of 7 X . On the other hand, again as in the ID case, 
the Qi terms are responsible for both lateral movement 
of the structure and shape asymmetry. These effects can 
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FIG. 6: Time evolution of relatively disordered ripples with 
mild wavelength coarsening (see also supplementary video).— 
Snapshots at increasing times: (a) t = 10; (b) t = 106; (c) 
t = 953 for Eq. J42J) with v = 0, = -0.1, v x = -1, v v = 

-o.i, = i, n y = 0.5, & = o.i, Ai x) = i, a' 1 ' = 5, a$ = 5, 

and ICij = 1. Top views (left column) and corresponding 
transverse cuts at y — L/2 (right column). Inset in (a) is its 
corresponding height autocorrelation. All units are arbitrary. 



be observed in Fig. [9] where we show the time evolution 
(as seen from a comoving reference frame) of transverse 
cuts of the surface for a given parameter condition. We 
have checked that, indeed, transforming back to a rest 
reference frame, the ripple velocity coincides, for times 
within linear regime, with that predicted by Eq. (|29[) . Al- 
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(c) 

FIG. 7: Time evolution of relatively ordered ripples with size- 
able wavelength coarsening (see also supplementary video).— 
Snapshots at increasing times: (a) t = 10; (b) t = 106; (c) 
t = 953 for Eq. (I42[l with the same parameters as in Fig. [6] 
except for = 0.1. Top views (left column) and corre- 
sponding transverse cuts at y = L/2 (right column). Inset in 
(a) is its corresponding height autocorrelation. All units are 
arbitrary. 




FIG. 8: Time evolution of relatively disordered ripples 
with sizeable wavelength coarsening (see also supplementary 
video).- - Ripple orientation as for typical large incidence an- 
gle conditions. Snapshots at increasing times: (a) t = 10; (b) 
t = 106; (c) t = 953 for Eq. J42) with v = 0,-y x = 0.1, v x = 1, 
v v = -0.95, a = -0.5, ii = 0.1, = 0.1, = 1.0, 
-^S = 0-5, A' 2 y = 5.0, and Kij = 1. Top views (left column) 
and corresponding transverse cuts at x — L/2 (right column). 
Inset in (a) is its corresponding height autocorrelation. All 
units are arbitrary. 



ready visual inspection of Fig. [5] suggests deviations from 
a uniform velocity for transverse ripple motion. This is 
a signature of nonlinear effects [specifically, due to rip- 
ple coarsening manifested by a non constant ripple wave- 
length l(t)], that are considered next. 



2. Non-linear regime 

For long enough times, non-linear terms have to be 
considered in order to understand the evolution of the 



16 



h(x) 




o 



50 



100 150 200 250 

X 



h(x)20 








50 



100 150 200 250 
X 



FIG. 9: In-plane non-uniform ripple motion as seen from the 
evolution of transverse cuts of the surface at y = L/2 for 
equally spaced times between t = and t = 1500. Results 
from the numerical integration of Eq. (|42[) with Vq = 0, v x = 



-1, u v = -0.1, 



A 



(2) 



50, A 



(2) 



Jx = £i = 0, 
= 5.0, and Ki,j 



2, Xx ' 1, Xy * — 5 

= 1. All units are arbitrary. 
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FIG. 10: Change in the direction of in-plane ripple motion 
as seen from the evolution of transverse cuts of the surface 
at y = L/2 for equally spaced times between t = and t = 
5000. Results from the numerical integration of Eq. (|42|l for 
parameters as in Fig. [9] except for £j = 4.5. All units are 
arbitrary. 



morphology. Those containing even derivatives are re- 
flection symmetric in x and, therefore, are not responsi- 
ble for lateral movement or any asymmetries of the pat- 
tern. On the other hand, we have checked that the terms 
corresponding to the coefficients £j indeed induce lateral 
motion of the pattern and asymmetry in the x axis. For 
the parameters considered in our simulations, positive 
values of £j induce a non-uniform lateral motion of the 
pattern towards positive x values. Since the contribu- 
tions of the £i nonlinearitics to the evolution of h increase 
in the non-linear regime, these can even induce a change 
in the direction of the pattern movement as observed in 
Fig. [ini where we plot the time evolution of a transverse 
cut of the surface. In this figure Sl^ = — 2 < 0, thus, as 
noted in the previous subsection, this induces a move- 
ment of the pattern towards negative x. These terms 
dominate during the linear regime but, as a result of the 
increase of the values of lower order surface derivatives, 
the £i = 4.5 > terms take over and change the direc- 
tion of lateral ripple motion towards positive x values. 
This example underscores the complex ripples dynamics 
induced by nonlinear effects, that should be taken into 
account in the discussion of the potential limitations of 
the current BH picture to quantitatively describe ripple 
motioni£i2i 

A simpler type of non-uniform ripple motion that has 
been reported experimentally corresponds to movement 
in a fixed direction, but with a non-uniform velocity, see 
e.g. 2lll65| . As mentioned above, this behavior corre- 



lates with the occurrence of wavelength coarsening (see 
below), and Eq. (|4"2]l is the first two-dimensional con- 
tinuum equation to describe it within the IBS context. 
As an example, in Fig. [TT] we show the (non-uniform) 



ripple velocity V(t) in the non- linear regime as a func- 
tion of time for the same simulations as shown in Fig. [9] 
Here the velocity is computed for a single surface mini- 
mum once the pattern is completely formed. At longer 
times the ripple velocity seems to reach a negligible value 
compatible with arrest of ripple motion. This might be 
related with a similar interruption of ripple coarsening 
that is illustrated below. 
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FIG. 11: (color online) Temporal evolution of the lateral rip- 
ple velocity, V(t), given by the numerical integration of Eq. 
(|42p in the non-linear regime for the same coefficients as in 
Fig. [9] A few representative error bars are given. All units 
are arbitrary. 



Non-linear terms containing derivatives that are 
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reflection-symmetric in x are responsible for the even- 
tual saturation of the ripple amplitude, and for the qual- 
ity and range of in-plane order of the ripple pattern. As 
checked in our simulations and described for the ID ef- 
fective equations studied in Sec. IV Bl and [4q . the larger 

the value of the ratio of \rj to A^ terms is, the longer 
is the coarsening process, and the more ordered the mor- 
phology becomes for the same total time (i.e., ion dose). 
This is shown in Figs. [T2l and 1131 in which the time evo- 
lution of the global surface roughness and of the lateral 
wavelength of the pattern are depicted for different val- 
ues of this ratio, denoted as r. In general, the roughness 
increases exponentially (linear instability regime), after 
which the nonlinearities are able to stabilize the system 
and induce slower growth for the roughness, W(t) finally 
reaching a time independent value. For very small r ra- 
tios, this stationary state seems to be reached earlier, and 
the intermediate slow (power-law) growth regime of the 
roughness is shorter. For larger values of r, this interme- 
diate regime has a wider duration, and can be more accu- 
rately described by a power law with the form W(t) ~ t@ 
for some effective value of the growth exponent (3. Note 
that, in the r — > oo limit (cquivalently, = 0), Eq. 
(|42p does not seem to have a stationary state, similarly 
to the conserved KS equation. 46 ' 66 Note, the growth ex- 
ponent for this case is^ 6 - (3 c ks = 1- The gradual chan ge of 
the duration of this intermediate power-law regime with 
physical parameters (that enter the value of the ratio r) 
and the different values for the effective growth exponent 
that can be obtained when trying to fit a power-law to 
such type of data, may account for the spread in the re- 
lated growth exponents experimentally reported in the 
context of ripple formation (see references in 

Regarding the quality and range of order in the rip- 
ple pattern at intermediate and long times, Fig. [12] al- 
ready shows that the morphology is more disordered (the 
roughness is larger) for smaller values of r. Moreover, for 
these cases, as seen in Fig. [13j the stationary value of 
the pattern wavelength is smaller, and is achieved ear- 
lier. A qualitatively similar behavior has been experi- 
mentally found in IBS of silicon targets under normal 
incidence conditions.— Note, the standard one-field con- 
tinuum equation (jAlj) corresponds to the r = limit, for 
which there is no coarsening and the system is roughest 
(the roughness being larger almost by an order of mag- 
nitude, as seen in Fig. [T2|) . Hence, such an equation was 
not able to account for the observed ripple coarsening and 
improved ordering, in marked contrast with the present 
Eq. (|42j) . As in the case of the roughness, for the opposite 
cKS-type limit r — > oo (A^ = 0), the ripple wavelength 
does not reach a stationary value. Rather, both the am- 
plitude and l(t) increase indefinitely, similarly to the cKS 
case for which Z(i) C KS ~ ^°' 5 i until a single ripple (with a 
parabolic cross-section) remains in a finite system.— 

Results obtained for the one dimensional anisotropic 
equation (|44p and for the ID and 2D isotropic 
counterparts^^ 6 - lead us to expect disorder to dominate 




FIG. 12: (color online) Temporal evolution of the global 
roughness, W(t), given by the the effective equation Eq. ((42 jt 
with do = 0, v x = — 1, v y — —0.1, Y-c = & = Qj = 0, Jdj = 1, 
A' 1 ' = 0.1, and A' 2 j = O.lr, for different values of r. The solid 

and dashed lines show the fit to power laws for A^ 1 ' = and 
r = 100 where W ~ i 1 ' 00 and W ~ t ' 71 , respectively. All 
units are arbitrary. 
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FIG. 13: (color online) Temporal evolution of the lateral 
wavelength of the pattern, l(t), given by the effective equation 
Eq. (|42|l for the same coefficients as Fig. 1121 The dashed line 
shows the fit to a power law for r = 100 where I ~ t°- 3S . The 
fit in the same region for A^ = yields I ~ t 0A8 (not shown). 
All units are arbitrary. 



the morphological features at the largest length and time 
scales in the system, as long as cancellation modes do 
not arise (namely, as long as A^ and Xy 1 * 1 have the same 
signs, see next section). Thus, we expect scale invari- 
ant morphologies and rough surfaces for much larger dis- 
tances than the pattern wavelength. The statistics of 
the surface fluctuations at these scales are expected to 
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be characterized by the critical exponents of some of the 
universality classes of kinetic roughening^ However, the 
case of the isotropic KS equation not being even com- 
pletely understood^ we can only conjecture, by analogy 
with the ID case, that the asymptotic scaling of Eq. ([44]) 
is in the 2D KPZ universality class. 

3. Cancellation Modes 

Eq. (|4"2"| can display cancellation modes (CM), analo- 
gously to its own ID counterpart, Eq. ([44]) . and to the 
anisotropic KS (aKS) equation ! 17 ' 22 Recall that CM in 
Eq. (|4"4"1) arise due to cancellation between the noncon- 
served (d x h) 2 and the conserved d 2 (d x h) 2 KPZ nonlin- 
caritics, and lead to (possibly) finite time blow-up of the 
solutions to the differential equation. We will refer to 
these as mixed KS (mKS) CM. In marked contrast, CM 
in the aKS system appear only when the coefficients of 
the two nonlinear terms A^ and Xy 1 ^ have different signs, 
and lead to a long time ripple pattern that is oriented 
along an oblique direction in the xy plane ; 17 ' 22 the sys- 
tem apparently supporting such type of solution for long 
times. We will denote these as aKS CM. 

Given the large parameter space of Eq. (|4"2"j) . the two 
types of CM mentioned can arise, and we consider sep- 
arately the conditions for appearance of each of them. 
Notice, it suffices to consider the nonlincarities that are 
reflection symmetric in x, as they are the only ones in- 
volved in the evolution (and putative blow up) of the 
ripple amplitude. 

a. mKS -type CM. The nonconserved and conserved 
KPZ nonlincarities in Eq. (|42j) read explicitly 

M[h] EE \^(d X hf + \ v l \dyh) 2 -\^dl(d X h) 2 (49) 

-XipUdyhf - \$d 2 y (d x hf - X$d 2 y (d y hf, 

whose Fourier transform reads 

T{M[h]) = Q x T[{d x h) 2 ] + Q x T[(d y h) 2 }, (50) 
where we have defined 

Qx = AW+A( 2 ifc 2 +A( 2 Jfc 2 , (51) 
Q y = AW+A^ + A^fc 2 . (52) 
Now, using (|4*3")) , we get 

Qi = Af } (l + ^fc 2 ) , (53) 

where we have assumed isotropy in the surface tension 
coefficients as done in Appendix [C| 722 = j2y = 72, and 
introduced A = <f>D/jo — 4>Reql2- As a function of system 
parameters, there are two possibilities: 

• If A > 0, then Q x ,Q y ^ 0, so that there are no can- 
cellations among nonconserved and conserved KPZ 
terms along any direction. This is the 2D general- 
ization of the analogous ID condition discussed in 
Section IV A 21 



• If A < 0, then cancellation occurs simultaneously 
in the x and y directions, for all Fourier modes 
on the circle |fc m _R-s| = (0/IAQ 1 / 2 , and we expect 
the solutions of Eq. (|4"2"|) to diverge for long times. 
However, as long as we are close to the instability 
threshold, the putative CM [being k m Ks ~ 0(1)] 
are outside the band of linearly unstable modes, so 
that no divergence occurs and Eq. (|4"2")) still provides 
a mathematically well-defined model. 

b. aKS- type CM. Even in the most favorable case 
(A > 0) considered in the previous discussion, there is 
still the possibility that cancellation takes place, not be- 
tween nonlincarities of different order (mKS type) but, 
rather, for specific directions on the xy plane, as in 
the aKS type. In order to assess such a possibility, 
we make the Ansatz^ 2 - that solutions are of the form 
h(x,y,t) = f(x — uy,t), and see how this reflects into 
the KPZ nonlinearities l|4"9"]) . Thus, 

N[h] = (AW + u 2 X^) jen 2 - ^j^Kf) 2 }"} , 

(54) 

where primes denote differentiation of / with respect to 
its first argument. As a consequence, exactly as in the 
aKS case, whenever the coefficients of the nonconserved 
KPZ nonlinearities, Xx and Xy , have different signs 
(as a result of their dependence on physical parameters), 
cancellation takes place for a Fourier mode that is ori- 
ented at an angle tan _1 (— X x ^ / Xy 1 ^) 1 / 2 with the x axis. 
Actually, for an appropriate choice of the function / ad- 
ditional cancellation may take place irrespective of the 
signs of the Kp coefficients, but such special cases are 
not generic. In Fig. [14] we show an example of the evo- 
lution of the morphology in case of cancellation modes 
oriented at 45° to the x axis. It is tempting to interpret 
the obliquely oriented ripples recently found^ on Si at 2 
keV as some type of cancellation mode of this aKS type. 

B. 2D dynamics: comparison to experiments 

Along the discussion on the detailed surface dynam- 
ics predicted by Eq. (|42]l . we have already pointed out 
relation to experimental features that are described by 
this new effective equation. In our discussion, we have 
assumed as a reference case that dependencies of coef- 
ficients a.jfcj on the physical parameters such as aver- 
age ion energy and flux, temperature, and characteris- 
tics of the distribution of energy deposition in the target, 
are as in BH-typc approaches for amorphizable targets. 
Within such an assumption, all dependencies of linear 
features on the latter are as in one-field models^ whose 
comparison to experiments has been reviewed in detail 
elsewhere Whenever discrepancies arise, some may 
be due to deviations of the actual collision cascade statis- 
tics from Sigmund's Gaussian formula, and this is a mat- 
ter of current active research. 37 ' 38 
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t = 50 

FIG. 14: Long time devel- 
opment of oblique ripple pat- 
terns at 45° to the x axis 
due to cancellation modes (see 
text). Top views of mor- 
phologies obtained by numer- 
ical integration of Eq. (|42|l 
with wo = 0, v x = -1, u y = 1, 
"fx = & = fii = 0, Ki,j — 
1, Ai 1 ' = 0.1, A' 1 ' = -0.1, 
Ag = 0.5, and A$ = -0.5 
at increasing times. All units 
are arbitrary. 



There are other features of our two-field model and 
of the ensuing Eq. (|42[) . that seem more robust to mod- 
ifications in the values of the parameters entering T eXl 
provided there is a morphological instability in the "sur- 
face tension" coefficients. Thus, the formation and fast 
stabilization of a stationary value for the thickness of 
the amorphous mobile layer has been assessed e.g. for Si 
both in Molecular Dynamics^ and in experiments, see [|| 
or @ (for energies of tens of keV) . Asymmetry in ripple 
cross sections has also been assessed both by microscopy 
(for Si, see [(|) or by techniques in reciprocal space {e.g. 
sapphire^) . Also wavelength coarsening has been pro- 
fusely documented, there being a large spread in the val- 
ues of the effective coarsening exponents, see references 
in Q and more recently e.g. '7Q\ for SiC. As for in-plane 
ripple motion, there is a smaller number of studies, al- 
though detailed studies (typically employing focused ion 
beams) are indeed available^ for Si and for glass^ the 
phenomenon having been reported also in atomistic sim- 
ulations of amorphous carbon targets^! 



VII. CONCLUSIONS 

In this paper we have considered in detail a two-field 
continuum description of nanopatterning by IBS in the 
the most general (anisotropic) case of oblique ion inci- 
dence. The explicit coupling of the dynamics for diffusing 
species at the surface with the evolution of the topog- 



raphy assumes exchange between mobile and immobile 
material at the upper boundary of the latter. This de- 
scription goes beyond the IBS case and can be employed 
as a phenomenological formulation of more general phe- 
nomena in Surface Science. In the particular case of IBS, 
this approximation leads to a stationary value for the 
density of diffusing surface species that is very quickly 
reached at, as compared with the dynamics of the sur- 
face morphology. This fact and the shape of the effective 
interface equation are robust to the specifics of the dis- 
tribution of energy deposition, while these reflect in the 
values of the model parameters. For the sake of refer- 
ence, in this work we have been assuming BH behavior 
for the latter. Under this assumption, we have explored 
the qualitative properties of the two-field model, and seen 
that it indeed provides a more comprehensive framework 
than previous continuum descriptions fulfilling most of 
the desiderata formulated in the Introduction. The adop- 
tion of a two-field description is not for mere mathemati- 
cal convenience; rather, it responds to its enlarged set of 
physical mechanisms (such as redeposition) entering the 
constitutive laws in a more natural way {e.g. the coupling 
between excavation and surface transport, or corrections 
to BH prediction for the linear velocity of transverse rip- 
ple motion). The parameter range for cancellation modes 
is partially restricted, and important experimentally ob- 
served nonlinear behavior, such as wavelength coarsening 
and non-uniform ripple motion, can be accounted for. 
Analogous conclusions can be reached at in the cases 
of normal incidence, or oblique incidence onto rotating 
targets, that arc considered in detail in a forthcoming 
worki22 

Still, some features that remain theoretically unex- 
plained, such as e.g. the lack of pattern formation for 
IBS of Si at near normal incidence in many experimental 
settings, may be due to inaccuracies in Sigmund's de- 
scription of the statistics of energy deposition through 
collisions inside the target^ Note that, once these are 
improved upon, the resulting effective parameters could 
be used in turn as inputs for the local excavation rate 
([3]). Thus, the two-field model is not restricted in princi- 
ple to Sigmund's Gaussian statistics. Moreover, through 
the mentioned natural coupling between transport and 
topography, and through the incorporation of redeposi- 
tion effects, it already provides (albeit, admittedly, in 
a simplified form that is susceptible to refinement) de- 
scription of additional material rearrangement due to ion 
impingement, that currently seems necessary for a the- 
oretical description of IBS with an improved predictive 
power i22i 
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APPENDIX C: RIPPLE WAVELENGTH AND 
ORIENTATION 



APPENDIX A: SIGMUND'S ENERGY 
DISTRIBUTION 

In [IH, the most general equation of motion for the 
local surface velocity reads 

d t h = -vo + 7xd x h + n x dlh + n 2 d x d 2 y h 

+ [& (d x h) + v x ] (d % x ti) + fc, (d x h) + v y ] (d 2 y h) 

-D xx d^h-D yy d^h-D xy d 2 x d 2 y h+^ {d x hf+^ (d y hf , 

(Al) 

where all coefficients are functions** of physical parame- 
ters such as 9, E, and the characteristics of Sigmund's 
Gaussian energy distribution, such as the average pen- 
etration depth, a, and the lateral widths of the distri- 
bution a and p,. In particular, the coefficients cor- 
respond to the so-called ion-induced effective smoothing 
terms whose terms have the same shape as those charac- 
teristic of surface diffusion, but are of a mere "geometric" 
origin related with describing the surface height at suffi- 
ciently high order in a Taylor expansion in height deriva- 
tives. For this (standard) choice, the values of coefficients 
a^i in the excavation rate T ex in ([3]) arc 

Q'O = V , ai x = — 7a,/Uo, U2x,y = ~Vx,y/Vo, 

a?,x = -Oi/u , ot Zy = -Q, 2 /v , auj = -Aj/t'o, 

CH5x,y = -£,x,y/vo, OtQ x ,y = -\ x ,y/v . (A2) 

APPENDIX B: FORMULAE FOR SEC. [ill] 
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In this appendix, we determine the ripple orientation 
and wavelength within linear theory. Our starting point 
is Eq. (|20[) which, neglecting 0(e 2 ) terms, can be rewrit- 
ten as 



Re(w k ) = -v x k 2 x - v y k 2 - K xx k% 



K~ h 4 — IC k 2 k 2 
'^xy^xvi 



-yy^y 



where 



-e(f)j R eq a 2x 
DR eq j2x 



+ eR eq j 

K-yy = DR eq j2y 



<\>Oi±xx 



v y = -e<jrj R eq a 2y 



V 70 



4>Reql2x )0L2 



eReqlO 



4>UAyy 



^ 7o 



(\>Reql2y ) OL 2 



JC xy = DR eq ("f2x + 72y) + tReqlO 
4>D 



2cf>a± 



■HI 



y n/ x n y> 

(CI) 



(C2) 
(C3) 



(C4) 



(C5) 



(a2x + a 2 y) + (j)Req(l2xO;2y + l2yOl2x) 

7o 



Note that the above 0{k A ) form is a consequence of our 
long wavelength approximation to o>k- However, the pre- 
cise shape of the coefficients is sensitive to the order con- 
sidered in the (independent) expansion in powers of e, 
resulting in Eqs. (|C2[) - (|C5[) . Thus, for instance, given 



that 72x and 72,, are positive, K xx , K. yy , JC xy are also al- 
ways positive to O(e ). However, the signs of their 0(e) 
contributions can change with experimental conditions, 
see e.g. [Hj|. 

The experimentally observed pattern is oriented along 
the direction which yields the maximum value of the real 
part of the dispersion relation, and its wavelength is asso- 
ciated to the wave vector, k = (k x , k y ), which maximizes 
Eq. (|C1[) . This vector must verify 



<9Re(u;k) 
dk x 



<9Re(wk) 
dk„ 



y) = 0. 



(C6) 



These conditions have the following independent solu- 
tions 



- ea 2j + e ^ a Atj k 2 k 2 



k = (0, 0) , ki 




OLlx 



(B3) 



AT 2 — 4AT K 

'^■xy ^'^xx'^yy 



K 2 
'^xy 



^K-xx^yy 



(C7) 



J2 a vtf 



J=x,y 



(B4) 



The solution k3 exists provided the arguments of the 
square roots of its two components are positive, and 
IC 2 V — AK, xx K yy ^ 0. Otherwise, the only solutions to 
(|C6p are k , ki, and k 2 . Moreover, as noted above, for 
large angles of incidence, v x is positive and ki is not 
defined. 
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Note, wave vector k3 implies a surface morphology 
with a periodicity that is aligned neither with the x nor 
with the y directions ("oblique ripples"). In order to 
study this solution, we write 



Re(wki) - Re(w k3 ) = 



xy"x 



Re(^ k2 ) - Re(a; k3 ) 



(J^xy^y 



(C8) 



vv I 

2 



4/Cyj, (JC^y 4:1C XX iCyy) 



(C9) 



so that the signs of (|C8|) , (|C9|) are given by that of 



(CIO) 



Moreover, straightforward algebra shows that the sign 
of the determinant of the Hessian matrix evaluated at 
k = k/3 is opposite to that of (|C10[) . In summary, 



(i) For positive values of (|C10|) . k3 is a saddle point, 
so that the absolute maximum of (|Clj) takes place 
either at ki or at k 2 . 



(ii) For negative values of (|C10|) . k3 provides the abso- 
lute maximum of (ICll). 



These results are valid for any value of e. To order O(e ), 
it is easy to see from (|C3[) - (|C5[) that K xy = K xx + K, yy , 
so that (|C10P equals {JC XX — JC yy ) 2 > and condition (i) 
above holds. This is the situation that occurs in most 
of the physical systems we will be considering, due to 
the smallness of the corresponding values of e. Higher 
order corrections are sensitive to high order details of the 
distribution of energy deposition. Thus, the sign of the 
0(e) term in (|C10|) is given, for isotropic thermal surface 
diffusion (j 2 x ="/2y), by the sign of 2a ixy - a ixx - a iyy . 
Hence, for specific choices of these effective smoothing 
coefficients it is conceivable that oblique ripples occur in 
our two field model for large e values, but we will not 
consider such situations in what follows. 

In order to decide which of the remaining solutions 
provide the absolute maximum of Re(wk), we finally sub- 
stitute the wave vectors given by (|C7|) into Eq. (|C1[) ; we 
obtain simply 



Re(wko) = 0, Re(w kl ) 



41C XX 



Re(w k2 ) 



4/C, 



Further discussion on the final orientation of the ripple 
structure can be found in the main text in Sec. IIIIB 21 



APPENDIX D: MULTIPLE-SCALE ANALYSIS 

In this appendix, we provide the details for the deriva- 
tion of the effective interface equation (|42jl . The setting 
is provided by formulae (|35|) to |38|) . For further conve- 
nience, we obtain a useful expression through addition of 
(El to (Efil) 



-(/)T ex +eD\7 2 R-e 3/2 d Tl R-e 2 dT,R. 

(Dl) 



e 3/2 d Tl h+e 2 d T2 h 



We will introduce the expansions (|39|) and (|40|) into 
(j4Tj) and (|D1|) and solve order by order in powers of e 1 / 2 . 

Order e° : To this order, as provided by Eq. (j4"Tj) , there 
is no contribution and we obtain 



Ro = 0. 



(D2) 



This means that the most important contribution to R 
vanishes near the instability threshold. Hence, as we al- 
ready noted when we obtained the planar solution, R will 
be only slightly altered from its planar value. 

Order e 1 / 2 : Again, there are no contributions to this 
order and from (1411) we obtain 



Ri = 0. (D3) 

Order e 1 : At this order, Eq. ([4Tj) reads 

i? 2 = -R eq y ■ feVAo), (D4) 

which yields the first correction to the expansion of R 
and depends on the curvatures of ho- As anticipated in 
the main text, R n contributions indeed depend of lower 
order h m terms. 

Order e 3 / 2 : From (1411) to this order we obtain 



Rs 



R eq V ■ (72V/11) + 4>R e qai x dxh 



(D5) 



We can substitute the previous values for the expansion 
of R into (|D1|) to finally obtain 



dnho = -c/rfoReqOtixdxho, 



(D6) 



which allows us to formally solve for ho. Note that this 
equation provides in-plane propagation as the leading 
contribution to h in the slow time scales, with the same 
velocity as predicted by the imaginary part of the linear 
dispersion relation, Eq. (I19p . 

Order e 2 . Following the previous scheme, from (|4ip 
we have 



i?4 



-R eq y ■ (l2Vh 2 ) + ^T ex (e 2 ) + — V 2 i? 2 , (D7) 



To 
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where we have used T ex (e 2 ) to denote the order e 2 con- 
tribution of T ex , given by 

r ea; (e 2 ) = joReq [ai x dxhi + V ■ (o^Vho) + V/io ■ (aeV/io 

(D8) 

From (|D1[) we obtain again a closed evolution equation 
for h to order e 2 , namely, 



d Tl hi + d T2 h = ~<f>T ex (e 2 ) + DV 2 R 2 . 



(D9) 



As we can see from the previous results, to obtain the 
temporal derivatives of h to order e"/ 2 we need to know 
R to order R n -2- Since we already know the value of the 
expansion of R up to R4 , we can obtain a closed evolution 
equation for h to order e 3 . 
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Order e 5 / 2 : From (|D1|) we have 
OtM + ^tM = -4>f ex (e 5 / 2 ) + DV 2 R 3 -d Tl R2, (D10) 

where we have again used Y ex {e > / 2 ) to denote the e 5 / 2 
order contribution to T ex . The time derivative of R2 can 
be computed by making use of Eqs. (|D4[) and (|D6|) . to 
get 

d Tl R 2 = <hoRl q aixV ■ ij2pxh ). (Dll) 

Order e 3 : Similarly to the previous step, from (jDip 
we now have 

^fta + 3t 2 /i 2 = e ,(e 3 ) + DV 2 R 4 - d Tl i? 3 - d T2 R 2 . 

JD12) 

Noting that r ex (e") do not depend explicitly on i?, and 
since we (formally) know the values of R2, R3, and R4 as 



functions of terms of the expansion of h, we can finally 
obtain a closed equation for the evolution of h up to 
order e 3 . To this end, using the relation between the 
slow temporal variables given by (|34p . we can write the 
time derivative of the expansion of h, Eq. (|4"D|) . as 



d t h = e 3 / 2 d Tl h + e 2 (d T2 h + d Tl hi) 
+ e 5/2 {&r 2 hi + d Tl h 2 ) + e 3 {d T2 h 2 + d Tl h 3 ) . (D13) 

Since all terms in (|D13p arc known, substituting Eqs. 
fD6jl . (|D9|l . l(DT0|) . and |D12|) into l(DT3|) . and using jDl| 
in order to simplify the C(e 3 ) contribution in (|D13|1 . we 
obtain 



d t h = e 3 / 2 {-07o^ eg ai,a^} + e 2 {-07o-R e9 [v ■ (arfh) + Vh • (aeV/i)] - £i? e? V 2 [V ■ (t^V/i)]} 



.5/2 



{-(jry R eq 



(4>DR eq l - (hoR 2 eq l2) Vd x h] } 



S X V ■ (a3_V/i) + (a x /i)V • (asV/i) 

{-07oi? eg V ■ (a4VV 2 ^) + V ■ [(4>DR eq I - <i> l0 R 2 eq ) 72J V [v ■ (arfh) + • (aeVft)] + 4>cfo R 2 eq a 2 lx d%h} , 

(D14) 



where I stands for the 2x2 identity matrix, we have 
again employed the expansion of h, Eq. (|40p . and we 
have neglected sixth order derivatives of h. 

In order to compare with the original model and the 
dispersion relation obtained in Sec. IIII1 we return to the 



original variables. Recalling that h = h' p + h, X = e 1 / 2 x, 
Y = e 1 / 2 ?/, and e = ao/^oReq), we finally obtain equa- 
tions (|42|) and (j43|) of the main text. 
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